//========================================================================================
// MerHashMap.hh
// ---------------------------------------------------------------------------------------
//
//----------------------------------------------------------------------------------------
// $Id: //CBT-depot/CBT++/include/MerHashMap.hh#12 $
//----------------------------------------------------------------------------------------
//
//========================================================================================
//.......1.........2.........3.........4.........5.........6.........7.........8.........9
//

#ifndef _MERHASHMAP
#define _MERHASHMAP

#include "CompBioTools.hh"
#include "Sequence3Bit.hh"
#include "Mer.hh"
#include "SimpleMerHash.hh"



//------------------------------------------
// * MerHashF
//------------------------------------------
// Computes the hash value of a mer.   This simple hash function computes the 
// XOR of bits of upper 16-mer with the bits of the lower 16-mer.  The XOR truth table has 
// 1/2 0s and 1/2 1s, so to 1st order (i.e. if bases were uniform) this hash function is 
// uniform. 
//
//
class MerHashF {
public:
  enum
    {
      bucket_size = 4,
      min_buckets = 8
    };

  inline size_t operator()(const Mer &M) const {
    register size_t  h = M.mUpper32 ^ M.mLower32;
    return(h);
  }
	
  bool equals_mer(Mer &a,Mer &b){
    if ((a.mUpper32 == b.mUpper32) && 
				(a.mLower32 == b.mLower32)) return(true);
    else return(false);
  }

	bool less_mer(const Mer &a,const Mer &b) {
		unsigned long long aw,bw;
		aw = a.ToWord();
		bw = b.ToWord();
		return(aw < bw);
	}

  //bool operator()(const Mer &a,const Mer &b) const{
	//	return(less_mer(a,b));
  //}

  bool operator()(const Mer &a,const Mer &b){
		return(less_mer(a,b));
  }
};

//------------------------------------------
// * MerSymmetricHashF
//------------------------------------------
// Hash function that does not distinguish between Mer and 
// Reverse-comp(Mer)
// 
class MerSymmetricHashF{
public:
  enum
    {
      bucket_size = 4,
      min_buckets = 8
    };

  inline size_t operator()(const Mer &M) const {
    Mer RCMer = M;
    RCMer.ReverseComplement();

    register size_t h;

    // Compute the hash for only the smaller mer, where
    // smaller is determined by treating Upper32/Lower32 
    // as an 8 byte long long. 

    // Clearly if upper32 is smaller then M is smaller
    if (M.mUpper32 < RCMer.mUpper32){
      h = M.mUpper32^M.mLower32;
    }else{ 
      // Maybe RC is clearly smaller than M
      if (RCMer.mUpper32 < M.mUpper32){
				h = RCMer.mUpper32 ^ RCMer.mLower32;
      }else{
				// The upper parts are equal so let the lower part 
				// decide. 
				if (M.mLower32 < RCMer.mLower32){
					h = M.mUpper32 ^M.mLower32;
				}else{
					h = RCMer.mLower32 ^ RCMer.mUpper32;
				}
      }
    }
    return(h);
  }

  bool operator()(const Mer a,const Mer b) const {
    // The mers match in the same sense. 
    if ((a.mUpper32 == b.mUpper32) && (a.mLower32 == b.mLower32)){
      return(true);
    }else{
      // see if the reverse complement matches. 
      Mer rcb = b;
      rcb.ReverseComplement();
      if ((a.mUpper32 == rcb.mUpper32) && (a.mLower32 == rcb.mLower32)){
				return(true);
      }else{
				return(false);
      }
    }
  }
};


//----------------------------------------------------------------------------------------
// * MerHashMapT
//----------------------------------------------------------------------------------------
// A template class to represent a MerHashMap.  A hash map stores a type of pair, so each 
// entry has a key and a value.  This is implemented as a template class to make it easy to
// substitute different hash functors (H) or different sattelite information (S).  
//
template<class H,class S>
class MerHashMapT: public hash_map<Mer,S,H>{
public:
};

// MerSet is built on the STL set class.  It's horribly 
// inefficient for the sorts of operations we want to do. 
// It's based on trees, and is efficient for set operations
// such as union and intersection, but for simple hashing it 
// is not a good choice. 
//----------------------------------------------------------------------------------------
// * MerSet
//----------------------------------------------------------------------------------------
// A set is a collection where the only thing you can do is Add, Remove, and ask whether
// the object is in the set.  This is distinct from a map, which maps a key to a value. 
// Having no sattelite information (values), this takes less space than Map if Set functionality 
// is all you need. 
//
// 
class MerSet: public hash_set<Mer,MerHashF>{
public:

  bool Contains(Mer key){
    MerSet::iterator it;
    it = find(key);
    if (it != end()) return(true);
    else return(false);
  }

  // Check if fwd or reverse complement of this sequence occurs.  
  bool ContainsFWDREV(Mer key);

  void ReadRepeatMers(istream &fin,int minCount);

  void Add(Mer key){
    insert(key);
  }

  void Add(string mer){
    Mer temp(mer);  // convert the string to a mer
    insert(temp);
  }
  void Add(Sequence3Bit &S);

  void Add(istream &fin);

  void Remove(Mer key){
    erase(key);
  }
};


//----------------------------------------------------------------------------------------
// * Mer2IntMap
//----------------------------------------------------------------------------------------
// A simple mer hash map that saves only the count of the mer. 
// 
// Obviously, we could also define an Int class and then instantiate 
// MerHashMap<Int> to get the same effect.   This is a common enough thing to do that 
// it's implemented here as it's own class, both as syntatic sugar and to save the 
// small extra cost of having an Int class to dereference.  
// 
//class Mer2IntMap: public hash_map<Mer,int,MerSymmetricHashF,eqMerSymmetric>{
// KJD POSSIBLE BUG ALERT: 
// KJD I *believe* that icc gets it's equal compare from the functor
// KJD and I *expect* that gcc will also.   Need to test whether this compiles
// KJD under gcc, and whether it works.  
class Mer2IntMap: public hash_map<Mer,int,MerSymmetricHashF>{
public:

  Mer2IntMap(){
    mTotalMersInHash = 0;
  }

  void Add(istream &inFile);
  void Add(istream &inFile,string Subset);
  void Add(istream &inFile,unsigned int Subset);

  void Add(Sequence3Bit &Seq);
  void Add(Sequence3Bit &Seq,string &Subset);
  void Add(Sequence3Bit &Seq,unsigned int Subset);

  void Add(Mer M);

  bool Contains(Mer M){
    Mer2IntMap::iterator it;
    it = (*this).find(M);
		  
    if (it != end()) return(true);
    else return(false);
  }

  bool BeginsWith(string &Subset,Sequence3Bit &Seq,int bIdx);

  int mTotalMersInHash;  
};

/*
//----------------------------------------------------------------------------------------
// * SymmetricMer2IntMap
//----------------------------------------------------------------------------------------
// Some sugar. 
class SymmetricMer2IntMap : public Mer2IntMapT<MerSymmetricHashF,eqMerSymmetric>{
public:
};

//----------------------------------------------------------------------------------------
// * Mer2IntMap
//----------------------------------------------------------------------------------------
// some sugar.
class Mer2IntMap : public Mer2IntMapT<MerHashF,eqMer>{
public:
};
*/


//----------------------------------------------------------------------------------------
// * BigMerLocation
//----------------------------------------------------------------------------------------
// Need another MerLocation for big sequences. 
// Ideally would like to be able to compare, say, 50Mbp pairwise. 
// 200Kbp vs. 20K (i.e. one BAC vs. all BACs).  
// 
// The most straightforward thing to do is to add a second int for 
// the BaseIdx.   
// 
// Listed from MSB to LSB: 
// 
// mData0:
// 1  bit   IsReverse.  Sense of the match 0 = fwd, 1 = reverse. 
// 19 bits  SourceRead (524,288 reads) 
// 12 bits  Unused.  (Just for now, to keep from having to rewrite SetSource, etc.) 
//
// mData1:
// 32 bits  BaseIdx
//
class BigMerLocation{
public:
  unsigned int mData;
  unsigned int mData1;
  
  BigMerLocation(){
    mData=0;
    mData1=0;
  }

  void InitFrom(BigMerLocation &BM){
    mData = BM.mData;
    mData1 = BM.mData1;
  }
  
  int SourceRead() const {
    unsigned int sread = mData;
    sread = sread & 0x7FFFFFFF; // Mask out the IsReverse bit. 
    sread = sread >> 12;  // Shift out the BaseIdx, leaving the read.      
    return((int)sread);
  }
  
  void SetSourceRead(int SourceRead){
    unsigned int s = (unsigned int) SourceRead;
    s = s << 12;                       // Shift it over 12 bits.   
    unsigned int mask = 0x80000fff;    // Mask to protect the top bit and the lower 12 bits. 
    mData = mData & mask;              // clear the old
    mData = mData | s;                 // Set it to s
  }
  
  int SourceBaseIdx() const {
    return(mData1);
  }
  
  void SetBaseIdx(int BaseIdx){
    mData1 = BaseIdx;
  }
  
  void SetIsReverse(){
    //unsigned int mask = 0x7fffffff;   // Clear the old high bit.  
    //mData = mData & mask;
    mData = mData | 0x80000000;   // Set high bit. 
  }
  
  bool IsReverse(){
    unsigned int mask = 0x8fffffff;
    bool isrev = mData & 0x80000000;
    return(isrev);
  }
};


//----------------------------------------------------------------------------------------
// * EmptyMerLocation
//----------------------------------------------------------------------------------------
// This is a mer location with no data part. 
// 
//
class EmptyMerLocation{
public:
  int SourceRead() const {
  }  
  void SetSourceRead(int SourceRead){
    // There is always only one source read when this class is used. 
  }  
  int SourceBaseIdx() const {
  }  
  void SetBaseIdx(int BaseIdx){
  }  
  void SetIsReverse()
  {    //unsigned int mask = 0x7fffffff;   // Clear the old high bit.  
  }  
  bool IsReverse(){
  }

  

};



//----------------------------------------------------------------------------------------
// * SingleSourceMerLocation
//----------------------------------------------------------------------------------------
// This is for big sequences where there is guaranteed to only be one source
// sequence.  For example, the dot plotting program MegaDot plots a single source
// against multiple queries.  For this program there is no need to save the index 
// of the source that is used.  By eliminating the SourceIdx, we can achieve greater 
// space efficiency.  
// 
// mData0:
// 1  bit   IsReverse.  Sense of the match 0 = fwd, 1 = reverse. 
// 31 bits  BaseIdx
//
//
class SingleSourceMerLocation{
public:
  unsigned int mData;
  
  SingleSourceMerLocation(){
    mData=0;
  }

  void InitFrom(SingleSourceMerLocation &BM){
    mData = BM.mData;
  }
  
  int SourceRead() const {
    return(1);
  }
  
  void SetSourceRead(int SourceRead){
    // There is always only one source read when this class is used. 
  }
  
  int SourceBaseIdx() const {
    unsigned int mask = 0x7fffffff;   // Clear the high bit.  
    unsigned int baseIdx = mData & mask;
    return(baseIdx);
  }
  
  void SetBaseIdx(int BaseIdx){
    unsigned int mask = 0x7fffffff;   
    BaseIdx = BaseIdx & mask; // Make sure high bit is clear. 
    mask = 0x80000000;

    mData = mData & mask; // Clear all but the high bit. 
    mData = mData | BaseIdx;
  }
  
  void SetIsReverse()
  {    //unsigned int mask = 0x7fffffff;   // Clear the old high bit.  
    //mData = mData & mask;
    mData = mData | 0x80000000;   // Set high bit. 
  }
  
  bool IsReverse(){
    unsigned int mask = 0x8fffffff;
    bool isrev = mData & 0x80000000;
    return(isrev);
  }
};




//----------------------------------------------------------------------------------------
// * MerLocation
//----------------------------------------------------------------------------------------
// A MerLocation is stored in the hash for each occurence of a mer in the sequence.  It 
// packs the information into a single int for space efficiency.  The SourceRead and the 
// BaseIdx of the Mer are encoded in an int as follows:
//
// Listed from MSB to LSB: 
// 
// 1  bit   IsReverse.  Sense of the match 0 = fwd, 1 = reverse. 
// 19 bits  SourceRead (524,288 reads) 
// 12 bits  BaseIdx    (4096 base positions)
// 
// Realistically, 16 bits, or 65K reads, is probably the most we can expect from a 1GB 
// machine, so 4 bits of the int could be used for other info if the need arises. 
//
class MerLocation{
public:
  unsigned int mData;

  MerLocation(){
    mData=0;
  }

  void InitFrom(MerLocation &M){
    mData = M.mData;
  }

  int SourceRead() const {
    unsigned int sread = mData;
    sread = sread & 0x7FFFFFFF; // Mask out the IsReverse bit. 
    sread = sread >> 12;  // Shift out the BaseIdx, leaving the read.  

    return((int)sread);
  }

  void SetSourceRead(int SourceRead){
    unsigned int s = (unsigned int) SourceRead;
    s = s << 12;                       // Shift it over 12 bits.   
    unsigned int mask = 0x80000fff;    // Mask to protect the top bit and the lower 12 bits. 
    mData = mData & mask;              // clear the old
    mData = mData | s;                 // Set it to s
  }

  int SourceBaseIdx() const {
    unsigned int mask = 0x00000fff; 
    unsigned int bidx = mData;
    bidx = bidx & mask; // Mask out everything but lowest 12 bits. 
    return((int)bidx);
  }

  void SetBaseIdx(int BaseIdx){
    unsigned int b = BaseIdx;
    unsigned int mask = 0xfffff000;
    mData = mData & mask;              // clear the old
    mData = mData | b;
  }

  void SetIsReverse(){
    //unsigned int mask = 0x7fffffff;   // Clear the old high bit.  
    //mData = mData & mask;
    mData = mData | 0x80000000;   // Set high bit. 
  }
  
  bool IsReverse(){
    unsigned int mask = 0x8fffffff;
    bool isrev = mData & 0x80000000;
    return(isrev);
  }
};


//----------------------------------------------------------------------------------------
// * MerLocationVecT
//----------------------------------------------------------------------------------------
// A vector of MerLocations.  The value part of MerHashMap.  Each mer in the hash has an 
// associated vector of MerLocations. 
//
template<class MerLoc>
class MerLocationVecT : public vector<MerLoc>{
public:

  void Print(){
    for(int i = 0;i< this.size();i++){
      int Source = (*this)[i].SourceRead();
      int Base = (*this)[i].SourceBaseIdx();
      bool Sense = (*this)[i].IsReverse();
      
      cout<<Sense<<"\t"<<Source<<"\t"<<Base<<endl;
    }
  }
};

//----------------------------------------------------------------------------------------
// * MerSearchResult
//----------------------------------------------------------------------------------------
// MerSearchResult stores a single result from a hash search.  
//
template<class MerLoc>
class MerSearchResult: public MerLoc{
public:

  MerSearchResult(){
    mQueryBaseIdx = 0;
    mMerHashFrequency = 0;
    mMerGlobalFrequency = 0;
  }

  unsigned int mQueryBaseIdx;
  unsigned int mMerHashFrequency;
  unsigned int mMerGlobalFrequency;
};


//----------------------------------------------------------------------------------------
// * MerSearchResultVec
//----------------------------------------------------------------------------------------
// MerSearchResultVec stores a list of MerSearchResults returned from a hash search. 
//
//
template<class MerLoc>
class MerSearchResultVec: public vector<MerSearchResult<MerLoc> >{
public:

  void Print(ostream &fout){
    fout<<"Source\tBase\tSense\tQueryBaseIdx\tmerFreq\tglobalMerFreq"<<endl;
    for(int i = 0;i< this.size();i++){
      int Source = (*this)[i].SourceRead();
      int Base = (*this)[i].SourceBaseIdx();
      bool Sense = (*this)[i].IsReverse();
      int merFreq = (*this)[i].mMerHashFrequency;
      int globalMerFreq = (*this)[i].mMerGlobalFrequency;
      fout<<Source<<"\t"<<Base<<"\t"<<Sense<<"\t"<<(*this)[i].mQueryBaseIdx;
      fout<<"\t"<<merFreq<<"\t"<<globalMerFreq<<endl;
    }
  }

	/*
	void ReadWheeler(istream &fin){
		MerSearchResult<MerLoc> temp;
		do{
			string line;
			GetLine(fin,line);

			int off = 0;
			int sourceIdx = NextInt(line,off);
			int queryIdx = NextInt(line,off);
			char sense = SplitNext(line,off);

			
			temp.SetSourceRead(0);
			temp.SetBaseIdx(baseIdx);
			if(){
				temp.SetIsReverse();
			}
			
			temp.mMerHashFrequency = 1;
			temp.mMerGlobalFrequency=1;

			push_back(temp);
		}while(!fin.eof());
	}
	*/

  void Print(){ Print(cout);};

  void SortBySourceBaseIdx();
  void SortBySourceRead();
  void SortByMerCount();

  void SortByDiagonal();
  void SortByRevDiagonal();

};


//class ReadMerHashMap: public MerHashMapT<MerHashF,MerLocationVec>{
//public:
//};

//class SequenceMerHashMap: public MerHashMapT<MerHashF,BigMerLocationVec>{
//public:
//};


//----------------------------------------------------------------------------------------
// * MerHashMapT
//----------------------------------------------------------------------------------------
// A template class to represent a MerHashMap.  A hash map stores a type of pair, so each 
// entry has a key and a value.  This is implemented as a template class to make it easy to
// substitute different hash functors (H) or different sattelite information (S).  
//
//template<class H,class S>
//class MerHashMapT: public hash_map<Mer,S,H>{
//public:

//};

//----------------------------------------------------------------------------------------
// * MerHashMap
//----------------------------------------------------------------------------------------
// The instantiated version of MerHashMapT, with MerHashF as the hash function and 
// MerLocationVec as the sattelite information.  
//
// This instantiated version also keeps a local copy of the sequences used to create the 
// hash.  For space efficiency, this local copy is a 3-bit encoding of the sequence 
// provided by the class Sequence3Bit (3 bit because we still want to know Ns and Xs). 
// 
template<class H,class S>
class MerHashMap : public hash_map<Mer,MerLocationVecT<S> ,H>{
public:

  int mMaxSeed;
  unsigned long  mKillFilteredCount;

  MerHashMap(){
    mKillFilteredCount = 0;
    mRepeatFilteredCount = 0;
    mMeanFilteredRepeatOccurence = 0.0;
    mMaxFilteredRepeatOccurrence = 0;
		mMaxMerLocations = 12;
  }


  MerHashMap(int bucketSize){
    mKillFilteredCount = 0;
    mRepeatFilteredCount = 0;
    mMeanFilteredRepeatOccurence = 0.0;
    mMaxFilteredRepeatOccurrence = 0;
    mBucketSize = bucketSize;
    MerHashMap<H,S>::resize(bucketSize);
		mMaxMerLocations = 12;
  }

  void SetMerSize(int msize){
    Mer::mMerSize = msize;
  }

  // Want to keep a local copy of the source sequences for reference. 
  // 3 bit encoding is used to keep this small, while still allowing 
  // N and X in input.  Note that storage is actually 10 bases per int32 w/ 2 bits unused.
  Sequence3BitVector  mSourceSeqs;

  int mBucketSize;

  int mRepeatFilteredCount;
  float mMeanFilteredRepeatOccurence;
  int  mMaxFilteredRepeatOccurrence;

  // Maximum locations allowed per mer. 
  int mMaxMerLocations;

  float mMeanMerFrequency;
  float mStDevMerFrequency;

  long   mTotalMersInHash;


  vector<string> mSourceNames;
  vector<int>    mSourceSizes;

  void ComputeMerStats();
  void RemoveRepeats(int repeatCutoff);
  

  void ScanToReadNumber(istream &fin,int ReadBegin,char * pBuf);
  

  void AddToHashMap(Sequence3Bit &Seq,int SeqIdx);
  void AddToHashMap(Sequence3Bit &Seq);

  void AddToHashMapSampled(Sequence3Bit &Seq,int SeqIdx, int skip);

  void AddToHashMap(Sequence3Bit &Seq,int SeqIdx,SimpleMerHash &SkipSet);
  void AddToHashMap(istream &inFile);  

  void AddToHashMap(istream &inFile,int ReadStart, int ReadEnd);


  void AddToHashMapSampled(istream &inFile,int ReadStart, int ReadEnd, int skip);
  void AddToHashMapThinned(istream &inFile,int ReadStart, int ReadEnd,
													 unsigned long long f, unsigned long long slice);

  void AddToHashMapThinned2(istream &inFile,
														unsigned long long f, unsigned long long slice);

  void AddToHashMapThinned(Sequence3Bit &Seq,int SeqIdx,unsigned long long f, unsigned long long slice);


  void AddToHashMap(istream &inFile,int ReadStart, int ReadEnd,SimpleMerHash &SkipSet);




  // Adds TempMer to hash if it is not already there. 
  // If it is already there, tacks on an entry for SeqIdx,BaseIdx 
  void AddLocation(Mer &TempMer,int SeqIdx,int BaseIdx);
  void AddMer(Mer &TempMer);

  void FindMerLocations(Sequence &s,int merStart,MerSearchResultVec<S> &Results);
  void FindRevMerLocations(Sequence &s,int merStart,MerSearchResultVec<S> &Results);

  void FindMerLocations(Sequence3Bit &s,int merStart,MerSearchResultVec<S> &Results,SimpleMerHash &M);
  void FindRevMerLocations(Sequence3Bit &s,int merStart,MerSearchResultVec<S> &Results,SimpleMerHash &M);

  void FindFwdRevMerLocations(Sequence &s,int merStart,MerSearchResult<S> &MLocations);

  void Print(ostream &fout){
    Sequence tempSeq;
    typename MerHashMap::iterator it;
    for(it = MerHashMap<H,S>::begin();it != MerHashMap<H,S>::end();it++){
      Mer tempMer = (*it).first;
      tempMer.ToSequence(tempSeq);
      fout<<tempSeq<<endl;
    }
  }

};


//========= BEGINING OF INCLUDE STUFF ===========================

#define SEQ_BUFFER_SIZE2 500000
extern char gseqbuf2[SEQ_BUFFER_SIZE2];


// KJD I'm a little confused by the template stuff.   It seems that 
// the only way to get the template functions here to work is to include this
// part of the file in MerHashMap.hh.   How painful.   
// 
// When I include this part of the file everything works OK. 
// When I compile this as part of a library, it seems to work, but 
// when I go to link the app the template functions fail to link. 
// 

// KJD: C++ general question needs to be answered: can you put template
// KJD: functions in a .cc and compile them?  I guess not, since it wouldn't
// KJD: know what classes to substitute in, and the compiler probably doesn't
// KJD: have a conditional compile feature. 


//------------------------------------------
// * ScanToReadNumber
//------------------------------------------
// Convenience method to scan through a stream, skipping reads until the specified read 
// number is found.   The first read in the file is read 0.  
// 
template<class H,class MerLoc>
void
MerHashMap<H,MerLoc>::ScanToReadNumber(istream &fin,int ReadBegin,char * pBuf)
{
  int readNum = -1;

  while((readNum <= ReadBegin) && (!fin.eof())){
    fin.getline(pBuf,SEQ_BUFFER_SIZE2);
    if (pBuf[0] == '>'){
      readNum++;		
    }
    if (readNum == ReadBegin) return;
  }		  
}


//------------------------------------
// * AddToHashMapThinned2
//------------------------------------
// Skips mers which are modulo f 
//
// Experimental version written explicitly for bangthin 
// Assumes a single very large source. 
// 
template<class H,class MerLoc>
void
MerHashMap<H,MerLoc>::AddToHashMapThinned2(istream &fin,
																					 unsigned long long f, unsigned long long slice)
{
  string line;
  // First line should be > 
  GetLine(fin,line);
  if (line[0] != '>'){
    cerr<<"ERROR AddToHashMapThinned2 stream does not begin with >"<<endl;
    exit(1);
  }else{
    mSourceNames.push_back(line.substr(1,line.size()-1));
  }

  Sequence S;
  long baseIdx = 0;
  do{

    GetLine(fin,line);

    if (S.size() > 0){
			
      // Extract the as-yet unused part of previous sequence. 
      int rstart = S.size()-Mer::mMerSize+1;
      int rend = S.size()-1;


      //	cout<<"Subsequence "<<rstart<<"\t"<<rend<<endl;
      S.GetSubSequence(S,rstart,rend);

      // Make a sequence out of the new line. 
      Sequence temp = line.c_str();
      // Concatenate it to the remainder from the previous
      S+= temp;

      //cout<<"S="<<S<<endl<<endl;

    }else{			
      S = line.c_str();

      //cout<<"S="<<S<<endl<<endl;
    }

    if (S.size() >= Mer::mMerSize){
      bool bContainsNX;
      int send = S.size()-Mer::mMerSize+1;
      //			cout<<"send = "<<send<<endl;
      for(int i = 0;i<send;i++){
				Mer tempMer(S,i,bContainsNX);			
				
				if (!bContainsNX){
					//	cout<<"Not NX"<<endl;
					// see if it's a mer we keep or a mer we ignore
					unsigned long long MerWord = tempMer.ToWord();
					//	cout<<"MerWord="<<MerWord<<endl;
					//cout<<"MerWord % "<<f<<" = "<<MerWord%f<<endl;
					if ((MerWord % f) == slice){	 
						// Adds TempMer to hash if it is not already there. 
						// If it is already there, tacks on an entry for SeqIdx,baseIdx
						//cout<<"Add "<<MerWord<<" at "<<baseIdx<<endl;
						(*this).AddLocation(tempMer,0,baseIdx);		
					}
				}
				baseIdx++;
      }

    }
  }while(!fin.eof());
}


//------------------------------------
// * AddToHashMapThinned 
//------------------------------------
// Skips mers which are modulo f 
//
template<class H,class MerLoc>
void
MerHashMap<H,MerLoc>::AddToHashMapThinned(Sequence3Bit &Seq,int SeqIdx,
																					unsigned long long f, unsigned long long slice)
{
  bool bContainsNX;
  for(int baseIdx = 0;baseIdx< (Seq.mBaseSize - Mer::mMerSize);baseIdx++){
    // Create a mer out of the next Mer::mMerSize bases starting at baseIdx
    // recording if N or X encountered. 
    Mer TempMer(Seq,baseIdx,bContainsNX);
		
    if (!bContainsNX){
      // see if it's a mer we keep or a mer we ignore
      unsigned long long MerWord = TempMer.ToWord();
			//			cout<<"slice= "<<slice<<"\tf= "<<f<<"\tMerWord= "<<MerWord<<" boolean: "<<((MerWord % f) == slice)<<endl;
      if ((MerWord % f) == slice){	 				
				// Adds TempMer to hash if it is not already there. 
				// If it is already there, tacks on an entry for SeqIdx,baseIdx
				(*this).AddLocation(TempMer,SeqIdx,baseIdx);		
      }
    }
  }
}


/*
  if (sizeof(long long) == sizeof(long)) {
  probe = mer % HASHSIZE;
  }
  else {
  // All these variables should be 32 bits except for mer
  long lo = mer & 0xFFFFFFFF;
  long hi = mer >> 32;
  probe = (lo ^ hi) % HASHSIZE;
  }

*/



//---------------------------------------
// * AddToHashMapThinned (fin,begin,end)
//---------------------------------------
// Experimental version that skips skip bases between hash add. 
//
template<class H,class MerLoc>
void
MerHashMap<H,MerLoc>::AddToHashMapThinned(istream &fin,int ReadBegin,int ReadEnd, 
																					unsigned long long f, unsigned long long slice)
{
  int reportNum = 5000;

  // returns the first line of 
  //  cout<<"\tScanning to first read in source..."<<flush;
  ScanToReadNumber(fin,ReadBegin,gseqbuf2);
  //cout<<"done.";;
  string temp = gseqbuf2;
  int newoffset;
  temp = GetNextWSDelimitedString(temp,0,newoffset);
  //cout<<"  First Read ="<<temp<<endl;
  int NumReads = ReadEnd - ReadBegin + 1;
    
  // Make space for new reads.
  // KJD Didn't clear mSourceSeqs because this function is "Add" 
  mSourceSeqs.resize(mSourceSeqs.size()+NumReads);

  int readNum = 0;
  int count = 0;
  do{
    //	 int seqIdx = readNum+ReadBegin;
    int seqIdx = readNum;
    
    // Probably want to make this optional. 
    //    mSourceSeqs[seqIdx].GetNextReadTrimLowQual(fin,gseqbuf2,SEQ_BUFFER_SIZE2);
    mSourceSeqs[seqIdx].GetNextReadTrimLowQual(fin,gseqbuf2,SEQ_BUFFER_SIZE2);

		//		cout<<"BUG: mSourceSeq.size="<<mSourceSeqs[seqIdx].mBaseSize<<endl;
    
    // KJD Might save doing the subtraction end()-begin() each time, with corresponding cost. 
    //    mSourceSizes.push_back(mSourceSeqs[seqIdx].size());

    // Want the number of bases in the sequences
    mSourceSizes.push_back(mSourceSeqs[seqIdx].mBaseSize);
    //	 cout<<dec<<seqIdx<<"\tAdd "<<mSourceSeqs[seqIdx].mBaseSize<<"\t"<<mSourceSeqs[seqIdx].size()<<endl;
  
    // This, on the other hand, is only referenced when saving overlaps, and could 
    // probably be replaced with a reference to mSourceSeqs. 
    mSourceNames.push_back(mSourceSeqs[seqIdx].mName);  
    
    count+=mSourceSeqs[seqIdx].mBaseSize;
    //	 cout<<"Add "<<mSourceSeqs[seqIdx].mName<<"\t"<<mSourceSeqs[seqIdx].mBaseSize<<endl;

    if (seqIdx == reportNum){
      int bcount = MerHashMap<H,MerLoc>::bucket_count();
      cout<<"\t** "<<seqIdx<<" reads ("<<count<<" bases) added.  Hash buckets: "<<bcount;
      cout<<"  HashSize ="<<MerHashMap<H,MerLoc>::size()<<" LoadFactor="<<((float)MerHashMap<H,MerLoc>::size()/(float)bcount);
      cout<<"\t"<<mSourceSeqs[seqIdx].mName<<endl;
      reportNum+=5000;
    }

    AddToHashMapThinned(mSourceSeqs[seqIdx],seqIdx,f,slice);	 
    readNum++;
  }while((!fin.eof()) && (readNum < NumReads));

}


//------------------------------------
// * AddToHashMapSampled
//------------------------------------
// Skips skip bases between samples. 
//
template<class H,class MerLoc>
void
MerHashMap<H,MerLoc>::AddToHashMapSampled(Sequence3Bit &Seq,int SeqIdx,int skip)
{
  bool bContainsNX;
  for(int baseIdx = 0;baseIdx< (Seq.mBaseSize - Mer::mMerSize);baseIdx+=skip){
    // Create a mer out of the next Mer::mMerSize bases starting at baseIdx
    // recording if N or X encountered. 
    Mer TempMer(Seq,baseIdx,bContainsNX);	 
    
    if (!bContainsNX){
      // Adds TempMer to hash if it is not already there. 
      // If it is already there, tacks on an entry for SeqIdx,baseIdx
      (*this).AddLocation(TempMer,SeqIdx,baseIdx);		
    }
  }
}


//------------------------------------
// * AddToHashMap (Seq)
//------------------------------------
// Adds a sequence, in 3-bit encoding, to the hash map, starting at SeqIdx.  That is, 
// the Mer::mMerSize bases starting at SeqIdx in Seq will be added to the hash.  Does nothing
// if the mer contains N or X. 
//
// KJD BUG:  Shouldn't this skip lower case also?   Does it??
//
template<class H,class MerLoc>
void
MerHashMap<H,MerLoc>::AddToHashMap(Sequence3Bit &Seq,int SeqIdx)
{
  bool bContainsNX;
  for(int baseIdx = 0;baseIdx< (Seq.mBaseSize - Mer::mMerSize);baseIdx++){
    // Create a mer out of the next Mer::mMerSize bases starting at baseIdx
    // recording if N or X encountered. 
    Mer TempMer(Seq,baseIdx,bContainsNX);	 
    
    if (!bContainsNX){
      // Adds TempMer to hash if it is not already there. 
      // If it is already there, tacks on an entry for SeqIdx,baseIdx
      (*this).AddLocation(TempMer,SeqIdx,baseIdx);		
    }
  }
}

//------------------------------------
// * AddToHashMap (Seq)
//------------------------------------
// Adds a sequence, in 3-bit encoding, to the hash map, starting at SeqIdx.  That is, 
// the Mer::mMerSize bases starting at SeqIdx in Seq will be added to the hash.  Does nothing
// if the mer contains N or X. 
//
// KJD BUG:  Shouldn't this skip lower case also?   Does it??
//
template<class H,class MerLoc>
void
MerHashMap<H,MerLoc>::AddToHashMap(Sequence3Bit &Seq)
{
  bool bContainsNX;
  for(int baseIdx = 0;baseIdx< (Seq.mBaseSize - Mer::mMerSize);baseIdx++){
    // Create a mer out of the next Mer::mMerSize bases starting at baseIdx
    // recording if N or X encountered. 
    Mer TempMer(Seq,baseIdx,bContainsNX);	 
    
    if (!bContainsNX){
      // Adds TempMer to hash if it is not already there. 
      // If it is already there, tacks on an entry for SeqIdx,baseIdx
      (*this).AddLocation(TempMer,0,baseIdx);		
    }
  }
}


//------------------------------------
// * AddToHashMap (Seq)
//------------------------------------
// Adds a sequence, in 3-bit encoding, to the hash map, starting at SeqIdx.  That is, 
// the Mer::mMerSize bases starting at SeqIdx in Seq will be added to the hash.  Does nothing
// if the mer contains N or X. 
//
// KJD BUG:  Shouldn't this skip lower case also?   Does it??
//
template<class H,class MerLoc>
void
MerHashMap<H,MerLoc>::AddToHashMap(Sequence3Bit &Seq,int SeqIdx,SimpleMerHash &SkipSet)
{
  int count=0;
  bool bContainsNX;
  for(int baseIdx = 0;baseIdx< (Seq.mBaseSize - Mer::mMerSize);baseIdx++){
    // Create a mer out of the next Mer::mMerSize bases starting at baseIdx
    // recording if N or X encountered. 
    Mer TempMer(Seq,baseIdx,bContainsNX);	 

    if (!bContainsNX){	 

      // SkipSet contains all of the repeat mers read in.  If this mer exceeds
      // MaxSeed, we omit it from the hash entirely and it can not be a seed 
      // for an overlap.  
      if (!SkipSet.Contains(TempMer,count)){
				(*this).AddLocation(TempMer,SeqIdx,baseIdx);		
      }else{
				// If it is in the skip set, see if it is below the max seed...
				if (count <= mMaxSeed){
					(*this).AddLocation(TempMer,SeqIdx,baseIdx);		
				}else{
					mKillFilteredCount++;
				}
      }
    }
  }
}

//---------------------------------------
// * AddToHashMap (fin)
//---------------------------------------
// Special version for the DotPlotter.  This will treat a multiple 
// sequence file as a single sequence.  Also, it will use 
// GetNextRead rather than GetNextReadTrimLowQual, so that it can 
// handle lower case letters in the input. 
//
template<class H,class MerLoc>
void
MerHashMap<H,MerLoc>::AddToHashMap(istream &fin)
{
  ScanToReadNumber(fin,0,gseqbuf2);

  string tempstr = gseqbuf2;
  int newoffset;
  tempstr = GetNextWSDelimitedString(tempstr,0,newoffset);
  
  // Adds all sequences as a single sequence. 
  mSourceSeqs.resize(1);
  
  int readNum = 0;
  int count = 0;
  Sequence3Bit temp;
  do{    
    temp.GetNextRead(fin,gseqbuf2,SEQ_BUFFER_SIZE2);    
    AddToHashMap(temp);
  }while(!fin.eof());
}

//---------------------------------------
// * AddToHashMapSampled (fin,begin,end)
//---------------------------------------
// Experimental version that skips skip bases between hash add. 
//
template<class H,class MerLoc>
void
MerHashMap<H,MerLoc>::AddToHashMapSampled(istream &fin,int ReadBegin,int ReadEnd, int skip)
{
  int reportNum = 5000;

  // returns the first line of 
  //  cout<<"\tScanning to first read in source..."<<flush;
  ScanToReadNumber(fin,ReadBegin,gseqbuf2);
  //cout<<"done.";;
  string temp = gseqbuf2;
  int newoffset;
  temp = GetNextWSDelimitedString(temp,0,newoffset);
  //cout<<"  First Read ="<<temp<<endl;
  int NumReads = ReadEnd - ReadBegin + 1;
    
  // Make space for new reads.
  // KJD Didn't clear mSourceSeqs because this function is "Add" 
  mSourceSeqs.resize(mSourceSeqs.size()+NumReads);

  int readNum = 0;
  int count = 0;
  do{
    //	 int seqIdx = readNum+ReadBegin;
    int seqIdx = readNum;
    
    // Probably want to make this optional. 
    //    mSourceSeqs[seqIdx].GetNextReadTrimLowQual(fin,gseqbuf2,SEQ_BUFFER_SIZE2);
    mSourceSeqs[seqIdx].GetNextRead(fin,gseqbuf2,SEQ_BUFFER_SIZE2);
    
    // KJD Might save doing the subtraction end()-begin() each time, with corresponding cost. 
    //    mSourceSizes.push_back(mSourceSeqs[seqIdx].size());

    // Want the number of bases in the sequences
    mSourceSizes.push_back(mSourceSeqs[seqIdx].mBaseSize);
    //	 cout<<dec<<seqIdx<<"\tAdd "<<mSourceSeqs[seqIdx].mBaseSize<<"\t"<<mSourceSeqs[seqIdx].size()<<endl;
  
    // This, on the other hand, is only referenced when saving overlaps, and could 
    // probably be replaced with a reference to mSourceSeqs. 
    mSourceNames.push_back(mSourceSeqs[seqIdx].mName);  
    
    count+=mSourceSeqs[seqIdx].mBaseSize;
    //	 cout<<"Add "<<mSourceSeqs[seqIdx].mName<<"\t"<<mSourceSeqs[seqIdx].mBaseSize<<endl;

    if (seqIdx == reportNum){
      int bcount = MerHashMap<H,MerLoc>::bucket_count();
      cout<<"\t"<<seqIdx<<" reads ("<<count<<" bases) added.  Hash buckets: "<<bcount;
      cout<<"  HashSize ="<<MerHashMap<H,MerLoc>::size()<<" LoadFactor="<<((float)MerHashMap<H,MerLoc>::size()/(float)bcount)<<endl;
      reportNum+=5000;
    }    

    AddToHashMapSampled(mSourceSeqs[seqIdx],seqIdx,skip);	 
    readNum++;
  }while((!fin.eof()) && (readNum < NumReads));

}


//---------------------------------------
// * AddToHashMap (fin,begin,end)
//---------------------------------------
// Adds the mer contents of the given filestream to the hashmap.  The file stream must be in 
// multiple-FASTA format.  This version adds only the reads from read number Begin to read 
// number End.  Mers with N or X or lower case are skipped. 
//
template<class H,class MerLoc>
void
MerHashMap<H,MerLoc>::AddToHashMap(istream &fin,int ReadBegin,int ReadEnd)
{
  int reportNum = 5000;

  // returns the first line of 
  //  cout<<"\tScanning to first read in source..."<<flush;
  ScanToReadNumber(fin,ReadBegin,gseqbuf2);
  //cout<<"done.";;
  string temp = gseqbuf2;
  int newoffset;
  temp = GetNextWSDelimitedString(temp,0,newoffset);
  //cout<<"  First Read ="<<temp<<endl;
  int NumReads = ReadEnd - ReadBegin + 1;
    
  // Make space for new reads.
  // KJD Didn't clear mSourceSeqs because this function is "Add" 
  mSourceSeqs.resize(mSourceSeqs.size()+NumReads);

  int readNum = 0;
  int count = 0;
  do{
    //	 int seqIdx = readNum+ReadBegin;
    int seqIdx = readNum;
    
    // Probably want to make this optional. 
    mSourceSeqs[seqIdx].GetNextReadTrimLowQual(fin,gseqbuf2,SEQ_BUFFER_SIZE2);
    
    // KJD Might save doing the subtraction end()-begin() each time, with corresponding cost. 
    //    mSourceSizes.push_back(mSourceSeqs[seqIdx].size());

    // Want the number of bases in the sequences
    mSourceSizes.push_back(mSourceSeqs[seqIdx].mBaseSize);
    cout<<seqIdx<<"\tAdd "<<mSourceSeqs[seqIdx].mBaseSize<<"\t"<<mSourceSeqs[seqIdx].size()<<endl;
  
    // This, on the other hand, is only referenced when saving overlaps, and could 
    // probably be replaced with a reference to mSourceSeqs. 
    mSourceNames.push_back(mSourceSeqs[seqIdx].mName);  
    
    count+=mSourceSeqs[seqIdx].mBaseSize;
    if (seqIdx == reportNum){
      int bcount = MerHashMap<H,MerLoc>::bucket_count();
      cout<<"\t"<<seqIdx<<" reads ("<<count<<" bases) added.  Hash buckets: "<<bcount;
      cout<<"  HashSize ="<<MerHashMap<H,MerLoc>::size()<<" LoadFactor="<<((float)MerHashMap<H,MerLoc>::size()/(float)bcount)<<endl;
      reportNum+=5000;
    }
    
    AddToHashMap(mSourceSeqs[seqIdx],seqIdx);	 
    readNum++;
  }while((!fin.eof()) && (readNum < NumReads));

}

//---------------------------------------
// * AddToHashMap (fin,begin,end)
//---------------------------------------
// Adds the mer contents of the given filestream to the hashmap.  The file stream must be in 
// multiple-FASTA format.  This version adds only the reads from read number Begin to read 
// number End.  Mers with N or X or lower case are skipped. 
//
template<class H,class MerLoc>
void
MerHashMap<H,MerLoc>::AddToHashMap(istream &fin,int ReadBegin,int ReadEnd,SimpleMerHash &SkipSet)
{
  int reportNum = 5000;

  // returns the first line of 
  cout<<"\tScanning to first read in source..."<<flush;
  ScanToReadNumber(fin,ReadBegin,gseqbuf2);
  cout<<"done.";;
  string temp = gseqbuf2;
  int newoffset;
  temp = GetNextWSDelimitedString(temp,0,newoffset);
  cout<<"  First Read ="<<temp<<endl;
  int NumReads = ReadEnd - ReadBegin + 1;
    
  // Make space for new reads.
  // KJD Didn't clear mSourceSeqs because this function is "Add" 
  mSourceSeqs.resize(mSourceSeqs.size()+NumReads);

  int readNum = 0;
  int count = 0;
  do{
    //	 int seqIdx = readNum+ReadBegin;
    int seqIdx = readNum;
    
    // Probably want to make this optional. 
    mSourceSeqs[seqIdx].GetNextReadTrimLowQual(fin,gseqbuf2,SEQ_BUFFER_SIZE2);
    
    // KJD Might save doing the subtraction end()-begin() each time, with corresponding cost. 
    //    mSourceSizes.push_back(mSourceSeqs[seqIdx].size());  

    mSourceSizes.push_back(mSourceSeqs[seqIdx].mBaseSize);  
    // This, on the other hand, is only referenced when saving overlaps, and could 
    // probably be replaced with a reference to mSourceSeqs. 
    mSourceNames.push_back(mSourceSeqs[seqIdx].mName);  
    
    count+=mSourceSeqs[seqIdx].mBaseSize;
    if (seqIdx == reportNum){
      int bcount = MerHashMap<H,MerLoc>::bucket_count();
      cout<<"\t"<<seqIdx<<" reads ("<<count<<" bases) added.  Hash buckets: "<<bcount;
      cout<<"  HashSize ="<<MerHashMap<H,MerLoc>::size()<<" LoadFactor="<<((float)MerHashMap<H,MerLoc>::size()/(float)bcount)<<endl;
      reportNum+=5000;
    }
    
    AddToHashMap(mSourceSeqs[seqIdx],seqIdx,SkipSet);	 
    readNum++;
  }while((!fin.eof()) && (readNum < NumReads));
}


//------------------------------------
// * AddLocation
//------------------------------------
// Adds TempMer to hash if it is not already there.  If it is already there, tacks on an 
// entry for SeqIdx,BaseIdx to the MerLocationVec for that Mer. 
//
template<class H,class MerLoc>
void 
MerHashMap<H,MerLoc>::AddLocation(Mer &TempMer,int SeqIdx,int BaseIdx){
  
	typename MerHashMap<H,MerLoc>::iterator it;

	
	//	cout<<"TempMer = ";TempMer.Print(cout);cout<<endl;

  it = (*this).find(TempMer);

  MerLoc MLocation;
  MLocation.SetSourceRead(SeqIdx);
  MLocation.SetBaseIdx(BaseIdx);

  //cout<<"BUG1:\t"<<MLocation.SourceRead()<<"\t"<<MLocation.SourceBaseIdx()<<endl;
  
  if (it != MerHashMap<H,MerLoc>::end()){
    
    // If it's already there, then just add another info to the list.
    // There is a hard limit to minimize effects of extreem outliers. 
    int merLocs = (*it).second.size();

		//cout<<"it != end() merLocs= "<<merLocs<<"\t < mMaxMerLocations = "<<mMaxMerLocations<<endl;

    if (merLocs < mMaxMerLocations){
      (*it).second.resize(merLocs+1);
      (*it).second[merLocs] = MLocation;
      //(*it).second.push_back(MLocation);
      mTotalMersInHash++;  // Used to compute a global average. 
    }
  }else{

		//		cout<<"it == end()"<<endl;

    // Otherwise, add a MerLocations structure to the list
    MerLocationVecT<MerLoc> MLocations;
    MLocations.push_back(MLocation);
    (*this)[TempMer] = MLocations;
  }
}

//------------------------------------
// * RemoveRepeats
//------------------------------------
// Removes the MerLocationVec entries for all mers that have more than repeatCutoff 
// separate occurrences.  
//
template<class H,class MerLoc>
void 
MerHashMap<H,MerLoc>::RemoveRepeats(int repeatCutoff)
{
  int TotalRepeatOccurences = 0;
  typename MerHashMap<H,MerLoc>::iterator it;
  for(it = MerHashMap<H,MerLoc>::begin();it != MerHashMap<H,MerLoc>::end();it++){
    // The value part of the hash is a MerLocationsVec structure, 
    // and it's size is the nubmer of distinct locations for this mer.
    int merCopies = (*it).second.size();
    if (merCopies > repeatCutoff){
      if (merCopies > mMaxFilteredRepeatOccurrence) mMaxFilteredRepeatOccurrence = merCopies;
      TotalRepeatOccurences += merCopies; 		
      (*it).second.clear();
      mRepeatFilteredCount++;
    }
  }
  
  mMeanFilteredRepeatOccurence = (float)TotalRepeatOccurences/(float)mRepeatFilteredCount;
}



//-------------------------------------
// * ComputeMerStats
//-------------------------------------
// Determine the mean and standard deviation of the mer frequency of what is actually in 
// the hash.  This may be different from what is in the sequence since there is a cutoff 
// on how many locations are recorded per position. 
//
template<class H,class MerLoc>
void 
MerHashMap<H,MerLoc>::ComputeMerStats()
{
  typename MerHashMap<H,MerLoc>::iterator it;
  double totalHits = 0;
  double numDistinctMers = 0;
  for(it = MerHashMap<H,MerLoc>::begin();it != MerHashMap<H,MerLoc>::end();it++){
    // The value part of the hash is a MerLocationsVec structure, 
    // and it's size is the nubmer of distinct locations for this mer.
    totalHits += (*it).second.size();
    numDistinctMers++;
  }
  mMeanMerFrequency = totalHits/numDistinctMers;
  double totalDeviation = 0;
  for(it = MerHashMap<H,MerLoc>::begin();it != MerHashMap<H,MerLoc>::end();it++){
    double diff = (*it).second.size() - mMeanMerFrequency;
    double deviation = diff * diff;
    totalDeviation +=deviation;
  }
  double variance = (1.0/(numDistinctMers - 1))*totalDeviation;
  mStDevMerFrequency = sqrt(variance);
}

//--------------------------------------
// * FindMerLocations
//--------------------------------------
// Finds all the locations of the mer specified by merStart in s. These locations are added 
// to the MerSearchResultVecVec.  
// 
template<class H,class MerLoc>
void
MerHashMap<H,MerLoc>::FindMerLocations(Sequence3Bit &s,int merStart,
																			 MerSearchResultVec<MerLoc> &Results,SimpleMerHash &KillSet)
{
  Mer tempMer(s,merStart); // 0.97
  typename MerHashMap<H,MerLoc>::iterator it;
  it = MerHashMap<H,MerLoc>::find(tempMer);  // 0.39 

  int count = 0;

  // Concatenate the new results onto the existing Results vector.
  if (it != MerHashMap<H,MerLoc>::end()) {
    MerLocationVecT<MerLoc> &foundLocations = (*it).second;
    register int lsize = foundLocations.size();

    // Repeat mers are marked by clearing the location vector part of the map. 
    // So, we only want to record hits that are to non-repeat samples.  
    if (lsize > 0){
      int oldSize = Results.size();	 
      Results.resize(oldSize+lsize);  // 0.67
      for(register int i = 0;i< lsize;i++){
				if (!KillSet.Contains(tempMer,count)){ 
					count = 0;
				}		  		  
				//		  Results[oldSize] = foundLocations[i];
				Results[oldSize].InitFrom(foundLocations[i]);
				Results[oldSize].mMerGlobalFrequency = count;
				Results[oldSize].mMerHashFrequency = lsize;
				Results[oldSize].mQueryBaseIdx = merStart;
				oldSize++;
      }
    }
  }
}

//----------------------------------------
// * FindRevMerLocations
//----------------------------------------
// Finds all the locations of the reverse complement of the mer specified by merStart in s 
// These locations are added to the MerSearchResultVec.  
// 
template<class H,class MerLoc>
void
MerHashMap<H,MerLoc>::FindRevMerLocations(Sequence3Bit &s,int merStart,
																					MerSearchResultVec<MerLoc> &Results,SimpleMerHash &KillSet)
{

  // KJD This is the origin of much later confusion.  The
  // way I sample and then RC is confusing to most people.  
  // It might be faster and less confusing to compute s.ReverseComplement()
  // and then just do a FindMerLocations with the Rev bit set. 
  // This will impact a lot of things, but is probably worth doing. 

  Mer tempMer(s,merStart);
  tempMer.ReverseComplement();

  typename MerHashMap<H,MerLoc>::iterator it;
  it = MerHashMap<H,MerLoc>::find(tempMer);

  int count = 0;

  if (it != MerHashMap<H,MerLoc>::end()) {
    // Is there a better way to concatenate vectors?  That's 
    // in essence what I want to do here. 
    MerLocationVecT<MerLoc> &foundLocations = (*it).second;
    register int lsize = foundLocations.size();

    // Repeat mers are marked by clearing the location vector part of the map. 
    // So, we only want to record hits that are to non-repeat samples.  
    if (lsize > 0){
      int oldSize = Results.size();
      Results.resize(oldSize+lsize);
      for(register int i = 0;i< lsize;i++){
				if (!KillSet.Contains(tempMer,count)){ 
					count = 0;
				}

				Results[oldSize].InitFrom(foundLocations[i]);
				Results[oldSize].mMerGlobalFrequency = count;
				Results[oldSize].mMerHashFrequency = lsize;
				Results[oldSize].mQueryBaseIdx = merStart;
				oldSize++;

				// KJD Why was I just copying mData?  Replaced with InitFrom() 2/28/03
				//		  Results[oldSize++].mData = foundLocations[i].mData;
      }
    }
  }
}

//-----------------------------------------
// * FindMerLocations
//-----------------------------------------
// Finds all the locations of the mer specified by merStart in s These locations are added 
// to the MerSearchResultVec.
// 
template<class H,class MerLoc>
void
MerHashMap<H,MerLoc>::FindMerLocations(Sequence &s,int merStart,
																			 MerSearchResultVec<MerLoc> &Results)
{
  Mer tempMer(s,merStart); // 0.97
  typename MerHashMap<H,MerLoc>::iterator it;
  it = MerHashMap<H,MerLoc>::find(tempMer);  // 0.39 

#if 0
  // KJD BUG
  
  for(int i = merStart;i< (merStart+32);i++){
    cout<<s[i];
  }
  cout<<endl;
  Sequence bugSeq;
  tempMer.ToSequence(bugSeq);
  cout<<"Mer: "<<bugSeq<<endl;
  // KJD BUG
#endif

  // Concatenate the new results onto the existing Results vector.
  if (it != MerHashMap<H,MerLoc>::end()) {
    MerLocationVecT<MerLoc> &foundLocations = (*it).second;
    register int lsize = foundLocations.size();

    // Repeat mers are marked by clearing the location vector part of the map. 
    // So, we only want to record hits that are to non-repeat samples.  
    if (lsize > 0){
      int oldSize = Results.size();	 
      Results.resize(oldSize+lsize);  // 0.67
      for(register int i = 0;i< lsize;i++){

				Results[oldSize].InitFrom(foundLocations[i]);

				//Results[oldSize] = foundLocations[i];

				Results[oldSize].mMerHashFrequency = lsize;
				Results[oldSize].mQueryBaseIdx = merStart;
				oldSize++;
      }
    }
  }
}

//--------------------------------------------
// * FindRevMerLocations
//--------------------------------------------
// Finds all the locations of the reverse complement of the mer specified by merStart in s. 
// These locations are added to the MerLocationVec.  
// 
template<class H,class MerLoc>
void
MerHashMap<H,MerLoc>::FindRevMerLocations(Sequence &s,int merStart,
																					MerSearchResultVec<MerLoc> &Results)
{
  Mer tempMer(s,merStart);
  tempMer.ReverseComplement();
  
  typename MerHashMap<H,MerLoc>::iterator it;
  it = MerHashMap<H,MerLoc>::find(tempMer);
  
#if 0
  // KJD BUG
  Sequence revs;
  for(int i = merStart;i< (merStart+32);i++){
    revs.push_back(s[i]);
  }
  revs.ReverseComplement();
  cout<<revs<<endl;
  Sequence bugSeq;
  tempMer.ToSequence(bugSeq);
  cout<<"RevMer: "<<bugSeq<<endl;
  // KJD BUG
#endif
  
  if (it != MerHashMap<H,MerLoc>::end()) {
    // Is there a better way to concatenate vectors?  That's 
    // in essence what I want to do here. 
    MerLocationVecT<MerLoc> &foundLocations = (*it).second;
    register int lsize = foundLocations.size();
    
    // Repeat mers are marked by clearing the location vector part of the map. 
    // So, we only want to record hits that are to non-repeat samples.  
    if (lsize > 0){
      int oldSize = Results.size();
      Results.resize(oldSize+lsize);
      for(register int i = 0;i< lsize;i++){
				Results[oldSize].InitFrom(foundLocations[i]);
				Results[oldSize].SetIsReverse();
				Results[oldSize].mMerHashFrequency = lsize;
				Results[oldSize].mQueryBaseIdx = merStart;
				oldSize++;
      }
    }
  }
}


//---------------------------------------------
// * cmpBySourceBaseIdx
//---------------------------------------------
// compare function for sort. 
//
template<class MerLoc>
bool cmpBySourceBaseIdx(const MerSearchResult<MerLoc> &a,const MerSearchResult<MerLoc> &b) {
  int aval = a.SourceBaseIdx();
  int bval = b.SourceBaseIdx();

  // break ties with the query base Idx. 
  if (aval == bval){
    return(a.mQueryBaseIdx < b.mQueryBaseIdx);
  }else{
    return(aval<bval);
  }
}


//---------------------------------------------
// * cmpByDiagonal
//---------------------------------------------
// Orders hits by which diagonal they are on.   
//
template<class MerLoc>
bool cmpByDiagonal(const MerSearchResult<MerLoc> &a,const MerSearchResult<MerLoc> &b) {
 
  int xa = a.SourceBaseIdx();
  int xb = b.SourceBaseIdx();

  int ya = a.mQueryBaseIdx;
  int yb = b.mQueryBaseIdx;

  int banda = ya - xa;
  int bandb = yb - xb;

  return(banda < bandb);
}


//---------------------------------------------
// * cmpByRevDiagonal
//---------------------------------------------
// Orders hits by which diagonal they are on.   
//
template<class MerLoc>
bool cmpByRevDiagonal(const MerSearchResult<MerLoc> &a,const MerSearchResult<MerLoc> &b) {
 
  int xa = a.SourceBaseIdx();
  int xb = b.SourceBaseIdx();

  int ya = a.mQueryBaseIdx;
  int yb = b.mQueryBaseIdx;

  int banda = xa-ya;
  int bandb = xb-yb;

  return(banda < bandb);
}


//----------------------------------------------
// * cmpBySourceRead
//----------------------------------------------
// compare function for sort. 
//
template<class MerLoc>
bool cmpBySourceRead(const MerSearchResult<MerLoc> &a,const MerSearchResult<MerLoc> &b) {
  int aval = a.SourceRead();
  int bval = b.SourceRead();
  return(aval<bval);
}


//----------------------------------------------
// * cmpBySourceRead
//----------------------------------------------
// compare function for sort. 
//
template<class MerLoc>
bool cmpByMerCount(const MerSearchResult<MerLoc> &a,const MerSearchResult<MerLoc> &b) {
  int aval = a.mMerGlobalFrequency;
  int bval = b.mMerGlobalFrequency;
  return(aval<bval);
}

//----------------------------------------------
// * SortByMerCount
//----------------------------------------------
// Sort the search results by base index. 
//
template<class MerLoc>
void 
MerSearchResultVec<MerLoc>::SortByMerCount(){
  stable_sort(this->begin(),this->end(),cmpByMerCount<MerLoc>);
}

//----------------------------------------------
// * SortBySourceBaseIdx
//----------------------------------------------
// Sort the search results by base index. 
//
template<class MerLoc>
void 
MerSearchResultVec<MerLoc>::SortBySourceBaseIdx(){
  stable_sort(this->begin(),this->end(),cmpBySourceBaseIdx<MerLoc>);
}


//----------------------------------------------
// * SortByDiagonal
//----------------------------------------------
// Sort the search results by the diagonal it is on. 
//
template<class MerLoc>
void 
MerSearchResultVec<MerLoc>::SortByDiagonal(){
  stable_sort(this->begin(),this->end(),cmpByDiagonal<MerLoc>);
}

//----------------------------------------------
// * SortByRevDiagonal
//----------------------------------------------
// Sort the search results by the diagonal it is on. 
//
template<class MerLoc>
void 
MerSearchResultVec<MerLoc>::SortByRevDiagonal(){
  stable_sort(this->begin(),this->end(),cmpByRevDiagonal<MerLoc>);
}

//----------------------------------------------
// * SortBySourceRead
//----------------------------------------------
// Sort the search results by source read. 
// 
template<class MerLoc>
void 
MerSearchResultVec<MerLoc>::SortBySourceRead(){
  stable_sort(this->begin(),this->end(),cmpBySourceRead<MerLoc>);
}



// =============== End of Include stuff

//#include <../src/MerHashMap2.cc>

#endif
